{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The roots of probabilistic graphical models go back to the 1980s, with a strong connection to Bayesian statistics. The story resembles that of neural networks: they have been around for over three decades and they need massive computational power. However, unlike in the case of deep learning, the requirements for computational resources remain out of reach. These models require sampling a distribution, and very often it is the Boltzmann distribution. Since quantum computers can give samples from this distribution, we can hope that quantum hardware can enable these models the same way graphics processing units enabled deep learning.\n",
    "\n",
    "# Probabilistic graphical models\n",
    "\n",
    "Probabilistic graphical models capture a compact representation of a joint probability distribution. For $\\{X_1,\\ldots,X_N\\}$ binary random variables, there are $2^N$ assignments. In a graphical model, complexity is dealt through graph theory. We get both an efficient treatment of uncertainty (probabilities) and of logical structure (independence constraints). The factorization of the probabilities happens along conditional independences among random variables. The definition is that $X$ is conditionally independent of $Y$ given $Z$ $(X\\perp Y|Z)$, if $P(X=x, Y=y|Z=z) = P(X=x|Z=z)P(Y=y|Z=z)$ for all $x\\in X,y\\in Y,z\\in Z$.\n",
    "\n",
    "The graph can be directed -- these are called Bayesian networks in general -- or undirected, in the case of Markov networks (also known as Markov random fields) [[1](#1)]. Graphical models are quintessentially generative: we explicitly model a probability distribution. Thus generating new samples is trivial and we can always introduce extra random variables to ensure certain properties. These models also take us a step closer to explainability, either by the use of the random variables directly for explanations (if your model is such) or by introducing explanatory random variables that correlate with the others.\n",
    "\n",
    "In a Markov random field, we can allow cycles in the graph and switch from local normalization (conditional probability distribution at each node) to global normalization of probabilities (i.e. a partition function). Examples include countless applications in computer vision, pattern recognition, artificial intelligence, but also Ising models that we have seen before: the factors are defined as degree-1 and degree-2 monomials of the random variables connected in the graph.\n",
    "\n",
    "The factorization is given as a sum $P(X_1, \\ldots, X_N) = \\frac{1}{Z}\\exp(-\\sum_k E[C_k])$, where $C_k$ are are cliques of the graph, and $E[.]$ is an energy defined over the cliques. If $P$ is a Boltzmann distribution over $G$, all local Markov properties will hold. The other way also holds if $P$ is a positive distribution.\n",
    "\n",
    "Let us define a Markov random field of binary variables. This will be a random Ising model over a three nodes. This will contain three cliques of a single node (the on-site fields) and two cliques of two nodes: the edges that connect the nodes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-19T20:10:30.684603Z",
     "start_time": "2018-11-19T20:10:30.190403Z"
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import dimod\n",
    "\n",
    "n_spins = 3\n",
    "h = {v: np.random.uniform(-2, 2) for v in range(n_spins)}\n",
    "J = {(0, 1): np.random.uniform(-1, 1),\n",
    "     (1, 2): np.random.uniform(-1, 1)}\n",
    "model = dimod.BinaryQuadraticModel(h, J, 0.0, dimod.SPIN)\n",
    "sampler = dimod.SimulatedAnnealingSampler()        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The probability distribution of a configuration $P(X_1, \\ldots, X_N) = \\frac{1}{Z}\\exp(-\\sum_i E[C_k])$ does not explicitly define the temperature, but it is implicitly there in the constants defining the Hamiltonian. So, for instance, we can scale it a temperature $T=1$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-19T20:10:32.696067Z",
     "start_time": "2018-11-19T20:10:30.687484Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD5CAYAAAAHtt/AAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAHupJREFUeJzt3Xd4nOWd7vHvMxoVq1jdNrZsyw0X3CQbGxIIGAihGQKYZrNZcpI9ORgMIaSycLIElhIgAWzMSc6mbWLTTe8Bh4QANpbk3rssV0lWLyNpnv1jRoPcpAE08065P9c11+iZGWluY123X37zFmOtRUREoofL6QAiIvL5qLhFRKKMiltEJMqouEVEooyKW0Qkyqi4RUSijIpbRCTKqLhFRKKMiltEJMq4Q/FD8/LybGFhYSh+tIhIzCopKam01ub39LqQFHdhYSErVqwIxY8WEYlZxphdwbxOoxIRkSij4hYRiTIqbhGRKKPiFhGJMipuEZEoo+IWEYkyISnuippm2ju8ofjRIiJxLyTFXd3o4dG/bgnFjxYRiXshG5UsWLqVpZsOhurHi4jErZDOuH/wzEr21jSH8i1EROJOSIrb7TIAHG5q4+bFpbRp3i0i0mtCUtxDclJJ8Jd36e4afvnWxlC8jYhIXApJcaclu/nRN0YH1v//Hzt4Z93+ULyViEjcCdmM+3+fOZxzx/QLrG9/bhW7q5pC9XYiInEjZMXtchkeuXoSg7L6AFDf0s5Ni0tpbe8I1VuKiMSFkO5VkpWaxILZRSQm+Obdaypq+c/XN4TyLUVEYl7ID3kvGpLNHReNDaz/++NdvLpqb6jfVkQkZoXlXCU3fKWQC8cPCKx/+sJqth9qCMdbi4jEnLAUtzGGB2dNZGhuKgCNng7mLiqlpU3zbhGRzytsZwfsm5LIE7OLSXL73nLj/np+/vK6cL29iEjMCOtpXccPyuQ/Zp4SWD+zopznS/aEM4KISNQL+/m4r5s2mG9OHhhY3/nSGjbtrw93DBGRqBX24jbG8J+XT2BEfhoALW1e5i4qobG1PdxRRESikiNXwElLdvPk9VNISfS9/bZDjdzx4hqstU7EERGJKo5duuzk/hnc+80JgfXLK/fy1PJyp+KIiEQNR685OWtKAddMHRxY/8er61hbUetgIhGRyOf4xYLvvuwUxgzIAMDT7uWmxaXUtbQ5nEpEJHI5XtwpiQksnFNMWlICALuqmvjJ86s17xYROQHHixtgeH46D1w5MbB+c+1+/vjRTucCiYhEsIgoboCZkwbyrdOHBtb3vbGBst2HHUwkIhKZIqa4Af794rFMGJQJQFuH5ebFZdQ0eRxOJSISWSKquJPdvnl3RoobgIqaZm5/dhVer+bdIiKdIqq4AQbnpPLIVZMC6/c2HuS3/9juYCIRkcgSccUNcP4pA/i3M4cF1g+9vYnlO6odTCQiEjkisrgBfnzBGIqHZAHQ4bXMe6qUyoZWh1OJiDgvYos7McHFgtnFZKcmAnCgrpXvP72SDs27RSTORWxxAwzM6sOvr5kcWH+4tZIF7291MJGIiPMiurgBzh7dj5tnjAysH31vM//cWulgIhERZ0V8cQN8/7xRnDY8BwBr4danyzhQ1+JwKhERZ0RFcbsTXDx+bRF56ckAVDZ4mPdUGe0dXoeTiYiEX1QUN0C/vik8ft1kXMa3Xr6jml+9u9nZUCIiDoia4gb4yog8bjvv5MB64d+2sXTjQQcTiYiEX1QVN8BNM0bytZPzA+vbnl1JRU2zg4lERMIr6orb5TL8+upJDOibAkBNUxs3Ly7F0655t4jEh6grboDc9GQWzC4iwT/wLttdw4NvbXQ4lYhIeERlcQNMLczhJxeMDqx/9+EO3lq738FEIiLhEbXFDfBvZw7nvLH9AusfPb+KXVWNDiYSEQm9qC5uYwyPXDWZguw+ANS3tHPT4lJa2jocTiYiEjpRXdwAmamJPDG7mMQE37x7bUUd976+3uFUIiKhE/XFDTBpcBZ3XjwusP7LJ7t5eWWFg4lEREInJoob4FunD+XiCScF1ncsWcO2Qw0OJhIRCY2YKW5jDA9cOYHC3FQAGj0d3LSolGaP5t0iEltiprgBMlISWThnCklu3x9r4/56fv7KWodTiYj0rpgqboBxA/vyi0tPCayfXbGH51aUO5hIRKR3xVxxA1xz6mCuKBoUWN/18lo27a93MJGISO+JyeI2xnDv5eMZ1S8dgJY2LzcuKqGxtd3hZCIiX15MFjdAapKbhXOK6ZOYAMD2Q43c8eIarNXFhkUkusVscQOM6p/BfVeMD6xfXrmXxct3O5hIROTLi+niBri8qIDrpg0OrO9+ZT1rK2odTCQi8uXEfHED/HzmKYw9qS8Ang4vcxeVUtfS5nAqEZEvJi6KOyUxgYVziklPdgOwu7qJHz+3WvNuEYlKcVHcAMPy0njwyomB9Vvr9vOHf+50LpCIyBcUN8UNcPHEk7jhK4WB9X1vbKB092HnAomIfAFxVdwAP7toDJMKMgFo91rmLS7jcKPH4VQiIsGLu+JOdiewYHYxfVN88+6KmmZ+8OxKvF7Nu0UkOsRdcQMMzknlkasnB9ZLNx3iN3/f7mAiEZHgxWVxA3x9XH++97XhgfXD72xi2fYqBxOJiAQnbosb4IffGM3UodkAdHgt854qo7Kh1eFUIiLdi+viTkxwMX92ETlpSQAcrG/l+0+vpEPzbhGJYHFd3AAnZfbh0WsmY3zXGubDrZXMf3+Ls6FERLoR98UN8LWT85k3Y2Rg/dh7W/hwS6WDiURETkzF7XfreSdz+vBcAKyFW58u40Bdi8OpRESOpeL2S3AZHrtuMvkZyQBUNXqYt7iM9g6vw8lERI6k4u6iX0YKj19bhMs/716+s5pH3t3sbCgRkaOouI9y+ohcbj9/dGD95N+28f7GAw4mEhE5kor7OG48awRnnZwfWN/2zCr2HG5yMJGIyGdU3Mfhchl+fc1kTspMAaC2uY2bF5fhade8W0Scp+I+gZy0JBbMLsLtH3ivLK/hgTc3OpxKRETF3a0pQ3P4yQVjAuvf/3MHb63d52AiEREVd4++e+Ywvj6uf2D9o+dWs6uq0cFEIhLvVNw9MMbw8KxJFGT3AaC+tZ25i0ppaetwOJmIxCsVdxAyUxNZOKeYpATff651e+u457X1DqcSkXil4g7SxIIs7rxkbGC9aNluXl5Z4WAiEYlXQRW3MeYFY8zFxpi4Lvp/OW0oF088KbD+2ZI1bD3Y4GAiEYlHwRbxk8BsYIsx5gFjzJieviEWGWN44IoJDMtLA6DJ08HcRSU0ezTvFpHwCaq4rbV/tdbOAYqBncC7xpiPjDHfNsYkhjJgpMlISeSJ2cUku33/6TYfaOCul9c6nEpE4knQow9jTC5wA/BdoAx4DF+RvxuSZBFs3MC+/OKyUwLr50v28OyKcgcTiUg8CXbGvQT4B5AKzLTWXmqtfcZaOw9ID2XASHX11MFcUTwosL7rpbVs2FfnYCIRiRfBbnH/l7V2nLX2fmvtPgBjTDKAtXZqyNJFMGMM935zPKP6+f7dam33ctOiUhpa2x1OJiKxLtjivvc4j33cm0GiUWqSmyevL6ZPYgIA2ysb+dmSNViriw2LSOh0W9zGmAHGmClAH2NMkTGm2H87G9/YJO6N7JfBfVeMD6xfXbWXvyzb7WAiEYl17h6e/wa+DyQLgF91ebweuCNEmaLO5UUFLN9xmKeW+wr7nlfXM7kgiwkFmQ4nE5FY1O0Wt7X2T9baGcAN1toZXW6XWmuXhCljVPj5zHGMO6kvAJ4OL3MXl1Db3OZwKhGJRT2NSq73f1lojPnB0bcw5IsaKYkJLJxTTHqy739iyqub+dFzqzTvFpFe19OHk2n++3Qg4zg36aIwL41fzpoYWL+z/gC/+3CHg4lEJBZ1O+O21v7Gf393eOJEv4smnMQNXynkjx/tBOCBNzdSNCSbKUOznQ0mIjGj2+I2xjze3fPW2lt6N05suOOisZSV17CqvIZ2r+XmxaW8fsuZ5KQlOR1NRGJAT6OSkh5uchxJbhcLrisis4/vNC77alv4wbMr8Xo17xaRL6+nUcmfwhUk1gzOSeVXV0/iO39aAcDfNh3iyQ+2cdOMkQ4nE5Fo19NeJY/67181xrxy9C08EaPXuWP7872zhgfWj7yziU+2VzmYSERiQU8H4PzZf/9wqIPEqh+eP5rSXYf5dOdhvBbmPVXGG7ecSX5GstPRRCRK9XQATon//gN85yY5DFQDH/sfkx4kJriYf11x4IPJQ/Wt3Pp0GR2ad4vIFxTsaV0vBrYBjwMLgK3GmAtDGSyWDMhM4dFrJmOMb/3Rtioee2+Ls6FEJGoFe3bAR4AZ1tqzrbVnATOAX4cuVuz52sn5zDtnVGA9//0t/H3zIQcTiUi0Cra4D1prt3ZZbwcOhiBPTLv13FF8ZUQuANbC959Zyf7aFodTiUi06WmvkiuMMVcA64wxbxhjbjDG/CvwKvBpWBLGkASX4bFriwIfTFY3epj3VCltHV6Hk4lINOlpi3um/5YCHADOAs4GDgE6hvsLyM9IZv51Rbj88+5Pdx7m4Xc2ORtKRKJKTwfgfDtcQeLJacNzuf380Tz0tq+wf/PBdk4dmsN54/o7nExEokGwe5WkGGNuMsYsNMb8vvMW6nCx7MazRnD26PzA+vbnVlFe3eRgIhGJFsF+OPlnYAC+K+J8gO+KOPWhChUPXC7Dr6+ezMDMFABqm9u4+akyPO2ad4tI94It7pHW2ruARv/5Sy4GJoQuVnzITkti/uxi3P6B96ryGu5/c4PDqUQk0gVb3J3X4KoxxowHMoHCkCSKM1OGZvPTC8cE1n/4507eXLPPwUQiEumCLe7fGmOygbuAV4D1wIMhSxVnvnPGMM7v8sHkj59fzc7KRgcTiUgkC6q4rbX/Za09bK39wFo73Frbr/PqOPLlGWN46KpJDM7pA0B9aztzF5XS0tbhcDIRiUTB7lWSa4yZb4wpNcaUGGMeNcbkhjpcPMnsk8jC2VNISvD9lazfV8cvXlvvcCoRiUTBjkqexneI+5XALKASeCZUoeLVhIJM7po5LrBevGw3L5VVOJhIRCJRsMWdY629x1q7w3+7F8gKZbB4df30IcycNDCwvuPFNWw9qD0vReQzwRb3UmPMtcYYl/92NfB6KIPFK2MM918xgeF5aQA0eTqYu6iUJk+7w8lEJFL0dJKpemNMHfA9YDHg8d+eBm4Lfbz4lJ7sZuH1xSS7fX89mw80cOdLa7FWF18QkZ6vgJNhre3rv3dZa93+m8ta2zdcIePRmAF9ueeb4wPrJaUVPLdij4OJRCRSBDsqwRhzqTHmYf/tklCGEp+rpw5m1pSCwPqul9eyYV+dg4lEJBIEuzvgA8Ct+A68WQ/c6n9MQuyey8Yzun8GAK3tXuYuKqW+pa2H7xKRWBbsFvdFwNettb+31v4euMD/mIRYn6QEnphTTGpSAgA7Khv52ZI1mneLxLGgRyUcuftfZm8HkRMb2S+d+6/47Jxer63ex18+2eVgIhFxUrDFfT9QZoz5ozHmT0AJcF/oYsnRLps8iDnThwTW97y2gdV7ahxMJCJO6bG4jTEG+BA4DVjiv51urX06xNnkKHddMo5TBvp25vF0+ObdtU2ad4vEmx6L2/qGqS9Za/dZa1+x1r5srd0fhmxylJTEBBbOKSYj2XfFuT2Hm/nh86s07xaJM8GOSj4xxpwa0iQSlKG5aTx01cTA+t31B/jdhzscTCQi4RZscc/AV97bjDGrjTFrjDGrQxlMTuyC8Sfxv746LLB+4M2NlOyqdjCRiIRTsMV9ITAcOAeYCVzivxeH/PTCMUwe7NvRp91ruXlxGdWNHodTiUg49HSukhRjzPeBH+Hbd7vCWrur8xaWhHJcSW4XT8wpJis1EYB9tS3c9sxKvF7Nu0ViXU9b3H8CpgJr8G11PxLyRBK0QVl9+NXVkwLrDzYf4skPtjmYSETCoafiHmetvd5/mbJZwJlhyCSfwzlj+nPj2SMC60fe2cTH26ocTCQiodZTcQd2ErbW6oTQEer2r5/MtMIcALwWbnm6jIP1LQ6nEpFQ6am4Jxlj6vy3emBi59f+83RLBHAnuJg/u4jctCQADtW3cutTK+nQvFskJvV0Pu4E//m4O8/J7e7ytc7HHUH6903hsWuLMMa3/nh7FY/9dbOzoUQkJD7PSaYkwp0xKo9bzx0VWM9fupUPNh9yMJGIhIKKO8bMO2cUZ4zMA8BauO2ZleyrbXY4lYj0JhV3jElwGR69djL9MpIBqG70MG9xGW0dXoeTiUhvUXHHoLz0ZOZfV0SCyzfwXrHrMA+/vcnhVCLSW1TcMWr68Fx+eP7owPo3f9/Ou+sPOJhIRHqLijuGfe9rwzlnTL/A+vZnV1Je3eRgIhHpDSruGOZyGR65ahKDsvoAUNfSzs2LS2lt73A4mYh8GSruGJedlsSC2UUkJvjm3av21HL/GxsdTiUiX4aKOw4UDcnmZxeODaz/+NFOXl+9z8FEIvJlqLjjxLe/WsgFpwwIrH/ywmp2VDY6mEhEvigVd5wwxvDLqyYyJCcVgIbWduYuKqWlTfNukWij4o4jfVMSWTinmCS37699w7467n51ncOpROTzUnHHmfGDMvn5zHGB9VPLy1lSusfBRCLyeam449DsaUO4bPLAwPrfX1zLlgP1DiYSkc9DxR2HjDHcd/kERuSnAdDc1sGNi0pp8uhaGSLRQMUdp9KS3SycM4WURN+vwNaDDdz54lqs1cUXRCKdijuOjR6Qwb3fnBBYLymr4JlPyx1MJCLBUHHHuVlTCrh6akFg/X9fWce6vbUOJhKRnqi4hbsvHc+YARkAeNq93LSolPqWth6+S0ScouIW+iQl8MScYtKSEgDYWdXET19Yo3m3SIRScQsAI/LTeeDKiYH162v28d8f73IwkYiciIpbAmZOGsi/nDY0sL739fWsKq9xMJGIHI+KW45w5yVjmTAoE4C2DsvcRaXUNmneLRJJVNxyhGR3Ak/MLiYjxQ1ARU0ztz+3UvNukQii4pZjDMlN5aFZkwLrv244yJ0vraWqodXBVCLSScUtx3XB+AF854xhgfWiZbs548Gl3PfGBg7Vq8BFnKTilhP6yQVjOHNUXmDd3NbBb/++nTMefJ+7X13HgboWB9OJxC8Titnl1KlT7YoVK3r950r4eb2Wt9ft57H3trBx/5FnEExyu7j21MH8n7NGMNB/QWIR+eKMMSXW2qk9vk7FLcHwei1/3XCAx9/fwtqKuiOeS0wwXDV1MDeeNYLB/ivsiMjnp+KWkLDWsnTTQR57b+sx+3i7XYYriwuYO2MEQ3PTHEooEr1U3BJS1lr+saWSx97bQsmuw0c8l+AyXDZ5IDfPGMnw/HSHEopEHxW3hIW1lo+3VfHYe1tYtqP6iOdcxnc05s0zRjKqf4ZDCUWih4pbwu6T7VXMf38L/9xadcTjxsBFE05i3jkjGTOgr0PpRCKfilscU7Krmsff28oHmw8d89w3TunPvHNGMd5/WL2IfEbFLY5bWV7D/Pe28N7Gg8c8d97Yfsw7ZxSTBmc5kEwkMqm4JWKs2VPL/Pe38M76A8c8d/bofOadM4opQ7MdSCYSWVTcEnHW761jwdItvLl2P0f/2p05Ko9554xi2rAcZ8KJRAAVt0SszQfqWfD+Vl5dvfeYAj9teA63nDuK04fnYoxxJqCIQ1TcEvG2Hmxg4dKtvLSyAu9Rv4anFmZzy7mjOGNkngpc4oaKW6LGzspGnli6lSVlFXQc1eCj+qUza0oBlxcNol/fFIcSioSHiluiTnl1Ewv/tpXnS/bQ1nHk76XLwNdOzmfWlALOG9uflMQEh1KKhI6KW6JWRU0z/+9v23ihdA9Nno5jnu+b4mbmpIHMmlLA5MFZGqVIzFBxS9RrbG3nzbX7eb6knE+2Vx/3NcPz05g1pYArigoYkKlRikQ3FbfElPLqJpaUVvBC6R52Vzcd87zLwFdH5jFrSgHfOGWARikSlVTcEpOstXy68zDPl5Tz+up9NB5nlJLsdlE8JJtpw3KYPjyHosHZ9ElSkUvkU3FLzGvytPP2uv08X7KHj7ZVHbNPeKfEBMPEgiymD8th2rAcpgzNJiMlMbxhRYKg4pa4UlHTzIule3ixrIJthxq7fa3LwPhBmUwr9BX5tGE5ZKUmhSmpyImpuCVu7a9tYfnOapbvqGLZ9mq2HGzo8XvGDMgIlPi0YTn0y9AHnRJ+Km4Rv6qGVj7dWc2yHdUs31HN+n11JxyrdBqelxaYkU8blssgXQxZwkDFLXICtc1tlOz6rMjX7Kml/ehj7o8yKKsP04fn+OfkuRTmpmr/cel1Km6RIDW2tlO2u4blO6r4ZEc1K8tr8LR7u/2efhnJvi1yf5GP6peOy6Uily9HxS3yBbW0dbCqvIblO6pZvrOaFTsP09x27G6HXWWnJnKq/8PO6cNyGTewLwkqcvmcVNwivaStw8vailpfkfvLvL6lvdvvyUh2M6UwO1DkEwZlkuR2hSmxRCsVt0iIdHgtG/fXsWz7Z0Ve3ejp9ntSEn0HBU0flsu0YTkUDcnS0Z1yDBW3SJhYa9l6sCHwYeeyHVUcqGvt9nsSEwyTCrICe61MGZpNerI7TIklUqm4RRxirWV3ddMRRV5e3dzt9yS4DOMH9vXvR57LtMIcMlN1dGe8UXGLRJC9Nc3+EvcdGNTT0Z3GwOj+GUwflsP04bmcWphDfkZymNKKU1TcIhHsUL3voKDOMt+4P4iDgvLTmD4sN3DOlYE6KCjmqLhFokhNk4cVOw+zfGc1y7ZXsXZv3TGXcTtaQXYfRvfPIDc9ibz0ZHLTk8kLfO27z05N0m6JUSTY4tanISIRICs1ifPG9ee8cf0BaGhtp2TXYZbvqGL5jmpWldfi6TjyoKA9h5vZc7j72bkxkJN6ZJl33uelJ5GblkxeRjK5ab7HdPrb6KDiFolA6cluzjo5n7NOzgd8BwX5ju6sZvnOKkp2HaalrfujOwGshapGD1WNHjjQ8/umJSUEttyP2IJPS/IXfDL5Gb7Cz+yTqKNFHaLiFokCKYkJnD4il9NH5AKj8LR72bCvjn21LVQ1tlLV4KGywXd/qKGVqoZWqho91DS1fa73afR00FjddNyrDB3N7TLkpJ2o4H33eWm+Lfzc9CSS3dqa7y0qbpEolOR2MWlwFpMGd/86T7uXw00eDtX7iryqofWogvdQ1dhKZb3vvq0j+M+82r2Wg/WtHKzvfp/1ThkpbvL9oxrfiOazUU1el38ActOT6Zvi1km8uqHiFolhSW4X/fum0L9vz+cXt9ZS19xO5RFb8K0cavAVfuCxRg+V9a3Ut3Z/2P/R6lvaqW9pZ3tl97tCAiQluAJb6r4teV/Rd27Bd87q89OTyU5LIjEhvk4noOIWEQCMMWSmJpKZmsiI/J5f39LWQXXj0SOaLlv1jZ9t6Vc3enrcS6YrT4eXfbUt7KttCer1WamJR4xqPtuCP+rD2PRk0pISon5rXsUtIl9ISmICA7P6BLU/uddrqWlu82/Bdy1434jmkH9U07lV33Sci0B3p6apjZqmNrYGldt11Igm+nanVHGLSMi5/B9k5qQlMap/Ro+vb/K0H/GB69Fb8F1n9dVNnh4PXuqqpc1LRU0zFTXd70oJvuuT5qT5ZvHH252ys/Bz05LIz0gO24nDVNwiEnFSk9yk5rgZnJPa42vbO7wcbmo74kPWIwv+s/tDDa09XiSjK6+FygYPlQ3B707ZuVdN56im61Z8b+1OqeIWkajmTnCRn5HsO5fLgO5fa62l0dNBZX3rcUc0vbI7ZVUTu6q+2O6UwVJxi0jcMMaQnuwmPdlNYV5aj6+PpN0pu1Jxi4icQCTtTtmViltEpBf0xu6U1zwY3HupuEVEHHC83SmvCfJ74+twIxGRGKDiFhGJMipuEZEoo+IWEYkyKm4RkSij4hYRiTIqbhGRKBOSq7wbYw4Bu3r9B4uIxLah1toeD98JSXGLiEjoaFQiIhJlVNwiIlFG5yqRqGKM6QDWdHnoaWvtA07lEXGCZtwSVYwxDdba9F7+mW5r7Rc/x6ZImGlUIjHBGLPTGHO3MabUGLPGGDPG/3iaMeb3xphPjTFlxpjL/I/fYIx5zhjzKvCOMcZljFlojFlnjHnNGPOGMWaWMeZcY8yLXd7n68aYJQ79MUUAFbdEnz7GmJVdbl3PhFlprS0GngR+6H/s34H3rbWnAjOAh4wxnZc+OR34V2vtOcAVQCEwAfiu/zmA94GxxpjOXbS+DfwhRH82kaBoxi3RptlaO/kEz3VuCZfgK2KA84FLjTGdRZ4CDPF//a61ttr/9RnAc9ZaL7DfGLMUwFprjTF/Bq43xvwBX6F/q/f+OCKfn4pbYknnxfs6+Ox32wBXWms3dX2hMWY60Nj1oW5+7h+AV4EWfOWuebg4SqMSiXVvA/OMMQbAGFN0gtd9CFzpn3X3B87ufMJauxfYC9wJ/DGkaUWCoC1uiTZ9jDEru6zfstb+tJvX3wM8Cqz2l/dO4JLjvO4F4FxgLbAZWAbUdnl+EZBvrV3/JbKL9ArtDijiZ4xJt9Y2GGNygeXAV621+/3PLQDKrLW/czSkCNriFunqNWNMFpAE3NOltEvwzcNvdzKcSCdtcYuIRBl9OCkiEmVU3CIiUUbFLSISZVTcIiJRRsUtIhJlVNwiIlHmfwAiVqNT2sipngAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "temperature = 1\n",
    "response = sampler.sample(model, beta_range=[1/temperature, 1/temperature], num_reads=100)\n",
    "energies = [solution.energy for solution in response.data()]\n",
    "fig, ax = plt.subplots()\n",
    "probabilities = np.exp(-np.array(sorted(energies))/temperature)\n",
    "Z = probabilities.sum()\n",
    "probabilities /= Z\n",
    "ax.plot(energies, probabilities, linewidth=3)\n",
    "ax.set_xlim(min(energies), max(energies))\n",
    "ax.set_xticks([])\n",
    "ax.set_yticks([])\n",
    "ax.set_xlabel('Energy')\n",
    "ax.set_ylabel('Probability')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, the conditional independences are already encapsulated by the model: for instances, spins 0 and 2 do not interact. In general, it is hard to learn the structure of a probabilistic graphical given a set of observed correlations in the sample $S$. We can only rely on heuristics. The typical way of doing it is to define a scoring function and do some heuristic global optimization. \n",
    "\n",
    "Once we identified or defined the graph structure $G$, we have to learn the probabilities in the graph. We again rely on our sample and its correlations, and use a maximum likelihood or a maximum a posteriori estimate of the corresponding parameters $\\theta_G$ with the likelihood $P(S|\\theta_G)$. This is again a hard problem.\n",
    "\n",
    "Applying the learned model means probabilistic inference to answer queries of the following types:\n",
    "\n",
    "-   Conditional probability: $P(Y|E=e)=\\frac{P(Y,e)}{P(e)}$.\n",
    "\n",
    "-   Maximum a posteriori:\n",
    "    $\\mathrm{argmax}_y P(y|e)=\\mathrm{argmax}_y \\sum_Z P(y, Z|e)$.\n",
    "\n",
    "This problem is in \\#P. Contrast this to deep learning: once the neural network is trained, running a prediction on it is relatively cheap. In the case of probabilistic graphical models, inference remains computationally demanding even after training the model. Instead of solving the inference problem directly, we use approximate inference with sampling, which is primarily done with Monte Carlo methods on a classical computer. These have their own problems of slow burn-in time and correlated samples, and this is exactly the step we can replace by sampling on a quantum computer.\n",
    "\n",
    "For instance, let us do a maximum a posteriori inference on our Ising model. We clamp the first spin to -1 and run simulated annealing for the rest of them to find the optimal configuration. We modify the simulated annealing routine in `dimod` to account for the clamping."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-19T20:10:32.992517Z",
     "start_time": "2018-11-19T20:10:32.705659Z"
    }
   },
   "outputs": [],
   "source": [
    "from dimod.reference.samplers.simulated_annealing import greedy_coloring\n",
    "\n",
    "clamped_spins = {0: -1}\n",
    "num_sweeps = 1000\n",
    "βs = [1.0 - 1.0*i / (num_sweeps - 1.) for i in range(num_sweeps)]\n",
    "\n",
    "# Set up the adjacency matrix.\n",
    "adj = {n: set() for n in h}\n",
    "for n0, n1 in J:\n",
    "    adj[n0].add(n1)\n",
    "    adj[n1].add(n0)\n",
    "# Use a vertex coloring for the graph and update the nodes by color\n",
    "__, colors = greedy_coloring(adj)\n",
    "\n",
    "spins = {v: np.random.choice((-1, 1)) if v not in clamped_spins else clamped_spins[v]\n",
    "         for v in h}\n",
    "for β in βs:\n",
    "    energy_diff_h = {v: -2 * spins[v] * h[v] for v in h}\n",
    "\n",
    "    # for each color, do updates\n",
    "    for color in colors:\n",
    "        nodes = colors[color]\n",
    "        energy_diff_J = {}\n",
    "        for v0 in nodes:\n",
    "            ediff = 0\n",
    "            for v1 in adj[v0]:\n",
    "                if (v0, v1) in J:\n",
    "                    ediff += spins[v0] * spins[v1] * J[(v0, v1)]\n",
    "                if (v1, v0) in J:\n",
    "                    ediff += spins[v0] * spins[v1] * J[(v1, v0)]\n",
    "\n",
    "            energy_diff_J[v0] = -2. * ediff\n",
    "        for v in filter(lambda x: x not in clamped_spins, nodes):\n",
    "            logp = np.log(np.random.uniform(0, 1))\n",
    "            if logp < -1. * β * (energy_diff_h[v] + energy_diff_J[v]):\n",
    "                spins[v] *= -1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Running this algorithm, we can obtain the most likely configuration:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-11-19T20:10:33.018780Z",
     "start_time": "2018-11-19T20:10:32.997312Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{0: -1, 1: 1, 2: -1}"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "spins"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Boltzmann machines\n",
    "\n",
    "A Boltzmann machine generates samples from a probability distributition $P(\\textbf{v})$ inferred from the data, where $\\textbf{v} \\in \\{0,1\\}^n$. The assumption is that this distribution lies on a latent space that can be paramerized by a set of hidden variables $\\textbf{h} \\in \\{0,1\\}^n$, such that $P(\\textbf{v})=\\sum_h P(\\textbf{v}|\\textbf{h})P(\\textbf{h})$. The joint probability distribution is modeled as a Gibbs distribution with the energy defined by an Ising Model: $P(\\textbf{v}, \\textbf{h})=\\frac{1}{Z} e^{-\\beta E(\\textbf{h},\\textbf{v})}$ and $E(\\textbf{h},\\textbf{v})=-\\sum_{i,j} W_{ij} h_i v_j$. It can then be shown that $p(\\textbf{h}|\\textbf{v})=\\sigma(W \\cdot \\textbf{v})$ and $p(\\textbf{v}|\\textbf{h})=\\sigma(W \\cdot \\textbf{h})$, where $\\sigma$ is the sigmoid function defined by $\\sigma(x)=\\frac{1}{1+e^{-x}}$.\n",
    "\n",
    "To train a Boltzmann machine, we look for the weights $W$ that maximizes the log-likelihood $L=\\sum_{\\textbf{v} \\in S} \\log(p(\\textbf{v}|W))$, where $S$ is the training set. This function can be optimized using regular gradient ascent: $W_{ij}^{(t+1)}=W_{ij}^{(t)} + \\eta \\frac{\\partial L}{\\partial W_{ij}}$. Computing the gradient $\\frac{\\partial L}{\\partial W_{ij}}$ is the hard part. Indeed, we can show that \n",
    "\n",
    "$$\\frac{\\partial L}{\\partial W_{ij}}=\\frac{1}{|S|} \\sum_{\\textbf{v} \\in S} \\mathbb{E}_{\\textbf{h} \\sim P(\\textbf{h}|\\textbf{v})}[h_i v_j] - \\mathbb{E}_{(\\textbf{h},\\textbf{v}) \\sim P(\\textbf{h},\\textbf{v})}[h_i v_j]$$.\n",
    "\n",
    "The first expectation value is easy to compute: it is equal to $\\sigma \\left( \\sum_j W_{ij} v_j \\right) v_j$. We only need to sum those expectation values over the dataset. This is called the positive phase, after its positive sign in the gradient.\n",
    "\n",
    "The second expectation value cannot be simplified as easily, since it is taken over all possible configuration $\\textbf{v}$ and $\\textbf{h}$. It would take an exponential amount of time to compute it exactly. We can use the exact same quantum sampling method as above to outsource this part of the calculation to a quantum processing unit and train Boltzmann machines."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# References\n",
    "\n",
    "[1] Koller, D., Friedman, N., Getoor, L., Taskar, B. (2007). [Graphical Models in a Nutshell](https://ai.stanford.edu/~koller/Papers/Koller+al:SRL07.pdf). In *Introduction to Statistical Relational Learning*, MIT Press. <a id='1'></a>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
